Phase-Field Model of Mode III Dynamic Fracture 



O 

o 

(N 



Alain Karma^, David A. Kessler^, and Herbert Levine'^ 
^Department of Physics, Northeastern University, Boston MA 02215 
^Department of Physics, Bar-Ban University, Ramat-Gan, Israel 
^Department of Physics, University of California, San Diego, La Jolla, CA 92093-0319 

(February 1, 2008) 

We introduce a phenomenological continuum model for mode III dynamic fracture that is based 
on the phase-field methodology used extensively to model interfacial pattern formation. We couple 
a scalar field, which distinguishes between "broken" and "unbroken" states of the system, to the 
displacement field in a way that consistently includes both macroscopic elasticity and a simple 
rotationaUy invariant short scale description of breaking. We report two-dimensional simulations 
that yield steady-state crack motion in a strip geometry above the Griffith threshold. 
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The nonequilibrium physics of dynamic fracture con- 
tinues to be a challenging topic of great interest . Re- 
cent efforts have been spurred by experimental findings 
regarding the breakdown of straight crack propagation 
(along with its associated smooth fracture surface) as 
the crack speed exceeds a threshold value. This insta- 
bility has been seen in a variety of materials, both crys- 
talline 1^ and amorphous (|,^, and it has been repro- 
duced in molecular dynamics simulations as well as 
with finite element schemes Q . 

For any material undergoing brittle fracture, linear 
continuum elasticity only provides an accurate descrip- 
tion of the displacements in regions that are not too close 
to the crack tip. The classic approach to this problem ||] 
has been to solve linear elastic equations, with bound- 
ary conditions providing the driving stresses, right up to 
this tip. This approach relies upon the assumption that 
the "process zone" inside of which the linear continuum 
equations break down is microscopic in size. The solu- 
tions have stress fields which become singular at the as- 
sumed tips, representing a finite flow of energy into the 
infinitesimally sized process zone. The velocity of the 
crack is then phenomenologically assumed to be given by 
some function of this energy flow rate. 

This approach has two main limitations from a physics 
perspective. Firstly, it does not provide insight into how 
the crack velocity is actually determined, e.g. how it 
depends on short-scale dissipation. Secondly, and more 
importantly, it fails to predict instabilities of the tip dy- 
namics. Thus, just as is the case in the well-studied prob- 
lem of dendritic solidification Q, one must supplement 
the macroscopic transport physics with a consistent, reg- 
ularizing microscopic theory on the tip scale in order to 
create a sensible theoretical framework. 

One method for accomplishing this task involves plac- 
ing the system on a lattice and allowing for the elas- 
tic forces to rapidly diminish at large atomic separation. 
Analytical and numerical studies of such 

models have shown that the details of the lattice struc- 
ture are critical for the tip dynamics. This is not sur- 
prising since in general the process zone scale is just the 
lattice spacing. Thus, these models are useful but can- 



not even qualitatively describe experiments in amorphous 
systems. What appears to be a more sensible approach 
for this class of systems is to construct a regularized con- 
tinuum model that maintains rotational symmetry even 
inside the process zone. Constructing such a theory is 
the aim of this paper Q . 

As a first step, we focus here on the simpler situation 
of mode III fracture for which the displacement u can be 
taken to be in a fixed direction (out of the plane) and 
hence can be represented by a scalar field u. Standard 
linear elasticity assigns the energy 



E = I dx - fie 
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with the strain e — Vu and an elastic constant fi. Allow- 
ing a material to fracture means that at large enough e^, 
the energy becomes strain independent, thereby elimi- 
nating the force. In an ideally brittle material, for exam- 
ple, this transition occurs immediately at some critical 
magnitude of the strain, Ec- Our basic idea involves rep- 
resenting the local state of the system, either unbroken 
with |e| < Ec or broken with |e| > ec, via a second "phase" 
field 4>{x,t). This field can be made to track the cor- 
rect state if it obeys a standard two-minimum Ginzburg- 
Landau equation with the relative energy of the two wells 
dependent on e^ ~ e^. Specifically, we choose 
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where yDiv((/>) = - 

ified later that has the properties g{Q) = 0, ; 
and ^'(0) = 5'(1) = 0. With these choices, the 
ima are always at (j) — and (j) — 1, with the 



and g is a function spec- 
g{l) = 1 
two min- 
absolute 



minimum shifting from 1 to as e ^ passes . 

To close the system, we need to specify how 
the elasticity equation. Note that the above 
follows from the relaxational dynamics, Tdt4> = 
where the energy E is now given by 



(j) affects 
equation 
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If we now interpret this E as the full potential energy 
including the elastic contribution, we see that the afore- 
mentioned properties of g{4>) will in fact eliminate the 
elastic force for large strain without having any effect at 
small strain where ~ 1 and thus g — 1- Consequently, 
our second equation is derived by varying this energy 
with respect to displacement, which yields 
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where we have allowed for both a Stokes drag term with 
coefficient b and a Kelvin viscosity rj. This equation com- 
pletes our model specification. 

Let us place this work in some perspective. Our ap- 
proach is similar in philosophy to, but very different in 
detail from, the work of Aranson and co-workers | |T^ ] 
who also derive a continuum two-field model for fracture. 
Most crucially, the physical interpretation of the phase 
field and hence the way in which it enters dynamically are 
completely different. As a result the present model avoids 
certain unphysical features of their model, e.g. the log- 
arithmic dependence of the crack opening on the system 
size. Within the traditional fracture community, several 
researchers have studied the effects of "softening" the 
elastic energy at large strain and compensating for the 
resultant instability in the equations by adding higher 
derivative terms. This approach, however, turns out to 
be difficult to extend to construct a continuum model 
where the fracture energy is independent of the system 
size and, at the same time, the strain is fully relieved in 
the bulk solid after passage of the crack. Moreover, it 
leads to a single fourth-order elasticity equation that is 
extremely hard to treat numerically. Finally, it is worth 
recalling that the original idea of representing differ- 
ent phases of a system via a field coupled to the macro- 
scopic dynamics, and thereafter using derivative terms 
in the phase field to regularize the problem, arose in the 
context of noncquilibrium crystal growth where it has 
become the method of choice [|l9|,^ for highly accurate 
computations of solidification microstructures. 

To understand how our model accounts for fracture, we 
start with the (one dimensional) snap-back of a stretched 
elastic band of size 2L after it breaks in the middle. Let 
us first consider the final time- independent cracked state; 
note that this is equivalently the asymptotic state for a 2- 
d crack once the tip has passed. This state is determined 
by solving the above equations with all time derivatives 
set to zero, with the boundary conditions u{±L) = ±A, 
4>{iL) = I. Note that 2 A is the total integrated strain 
that is conserved by the dynamics. Moreover, both e{y) 
and (j){y) are symmetrical about y = 0; thus we only need 
to find a solution in the interval < y < L. The elas- 
ticity equation, Eq. ^, requires e(y) = dyU — eo/g{(/){y)). 
Substituting this form into Eq. yields 



= D^c^" - V{jw{4>) 
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FIG. 1. Plots of the effective potential for 1-d static crack 
profiles — 1 and — 1/2). 



From now on, we rescale lengths to make = 1. Eq. 

can be thought of as the equation of motion of a ball 
rolling in an effective potential 
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Veff[4>) = -Vdw + 



(6) 



A schematic picture of this potential is shown in Fig. 1. 
The solution of interest corresponds to rolling in a "time" 
L from the top of the hill at = 1 to the turning point 
4>* located near = 0; this turning point exists because 
e^/g{4>) becomes large and positive for small 4>. 
The asymptotic steady-state crack is thus given by 
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The two unknown constants eq and Eq are fixed by the 
requirements that the above equation yields (/) = 1 at 
y ~ L and by the overall integrated strain constraint 
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For this solution to be physically acceptable, almost all 
of the displacement must occur in the crack, thereby re- 
lieving the strain in the bulk, which imposes a constraint 
on the form of the function g. To see why, consider the 
large L limit (where with = 1 our length unit is the 
process zone scale) and note that the standard Griffith 
criterion for fracture necessitates the scaling A ^ VX. 
There is a contribution to the integral on the left hand 
side of Eq. ^ that arises from (p of order (f>*, which is 
close to zero. If we choose g to vanish near (p — a,s 
a power law g ~ 0^^", it is clear from the form of the 

2 

effective potential that cj)* ~ eg+" . Given this, the local 

contribution to the integral scales like eg . This will 
match the right hand side with the choice 
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FIG. 2. Results of Id crack computations, (a) The strain 
at the crack center vs. time (b) asymptotic profiles of the 
strain and phase fields. L — 5, eo = 0.045 

eo - (9) 

Hence, as long as a is positive, eo will go to zero at 
large L fast enough such that the local contribution, i.e. 
that of the crack, to the overall displacement is dominant 
compared to the bulk contribution which scales as eoL. 
We note that the least residual stress occurs in the limit 
a — > 0+, which gives rise to an exponential decrease as 
a function of system size. Finally, the fracture energy 7 
remains finite as L gets large, as is also required for a 
sensible theory. In this limit, eo ^ and Eq fj,e^/2, 
and it is easy to derive the expression 

7 = %/2 / d4>Jfiel/2-VEFF{4>; £0 = 0) (10) 
Jo 

We now turn to the Id time-dependent problem. We 
choose g = 40'^ — so that a = 1 and eo should scale 
as L~'^/^; also we pick Cc = 1/2, /i = 1. Fig. 2 shows 
the time development of the strain at the origin, for a 
completely overdamped system, p = 0, & = 1, = 0. 
We see that, starting from an initial large strain region 
near the origin, the system proceeds to fully crack and 
approach the aforementioned asymptotic state. Notice 
that for this case of overdamped dynamics, there are two 
time scales visible in the strain relaxation; a fast time 
scale during which the crack (say in the (j) field) develops 



and a slower one during which strain is drained from 
the bulk into the crack. For a typical underdamped case, 
(p = 1, & = 0, r/ = 0.2), the second-stage relaxation to the 
steady-value involves, as it must, damped oscillations. In 
either case, in the second stage, the (j> field is slaved to 
the displacement making the system of equations very 
stiff. We coped with this difficulty using an implicit time- 
stepping scheme. 

The real test of any fracture model comes in two (or 
higher) dimensions. Now, the crack tip must advance 
by providing enough stress to strain new material be- 
yond a critical extension. We have carried out a pre- 
liminary simulation study of our model in a 2-d strip 
geometry (with the edges of the strip at y = ±L) using 
a standard Crank-Nicholson alterning-direction-implicit 
scheme . We used the initial condition corresponding 
to a strained solid with (t>{x, y) = I and u{x, y) = Ay/ L, 
which must produce a crack propagating along the x- 
direction above the standard Griffith threshold, A > Ac, 
where here Ac ~ a/StL and 7 = 0.3808... is given by Eq. 
[To| for the present model parameters (/z = 1, ec = 1/2). 
A typical time sequence for a stably propagating crack 
is presented in Fig. 3. We have checked that our results 
are reasonably independent of the discretization scale dx, 
such that we are truly seeing the results of the continuum 
regularization of the tip scale dynamics. 

In Fig. 4, we present the steady-state crack velocity as 
a function of the driving. The crack propagates above a 
critical value of the drive that is within a few percent of 
the analytically predicted value Ac = -/27L for the large 
L limit even though L is not so large (L = 10) in the sim- 
ulations. It is important to recognize that 1^(A) cannot 
be obtained from the usual continuum theory without 
additional assumptions; here it follows directly from the 
fact that we have a consistent theory at both the macro- 
scopic and microscopic scales. One can obtain similar re- 
sults from lattice models of fracture [p| jic|j2^ , at the price 
of introducing lattice scale instabilities [ ^2yi2| . These in- 
stabilities are connected to spatial period-doubling, when 
the times for breaking alternating diagonal bonds in a 
hexagonal lattice become unequal. Thus, they have no 
direct relevance for amorphous systems without any un- 
derlying crystallinity. Here, at least at the drivings we 
have investigated so far, we find no such instability, and 
steady cracking appears to persist to rather high dis- 
placements. In terms of experiment, this seems to imply 
that the observed transition at moderate displacement 
to more complex tip dynamics is presumably connected 
to the mode I nature of the fracture geometry, with its 
propensity to branch in the direction of the off-axis stress 
maximum as suggested originally by Yoffc ||2^. An ex- 
tension of the present approach to mode I is presently 
underway to test this hypothesis. 

In terms of physics, our approach leads to the introduc- 
tion of a new time scale, t, connected to the relaxational 
rate of the phase-field. At any non-zero r, this relaxation 
is a possible source of tip-scale dissipation and hence it 
can affect the crack propagation. Decreasing this pa- 
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FIG. 3. 2-d simulation snapshots; pictured here is in 
grey-scale (0=black, l=white) over 1/2 the computational 
domain of size 20x20; here p = 1, b = 0, rj = .2, A = 2.81, 
timestep=.02 and grid spacing = .05 
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FIG. 4. Velocity vs. driving from 2-d simulations. 

rameter in the simulations indeed increases the velocity, 
consistent with convergence to a finite limit as t — > 
with 0(t) corrections. This indicates that the velocity 
is predominantly limited here by the rate at which the 
strain is drained from the bulk into the crack. 

An important numerical issue for future consideration 
concerns the sharpness of the strain profile. To fully re- 
solve the spatial scales in the 2-d crack is a daunting task 
which probably cannot be accomplished by sticking with 
a fixed computational grid. A related issue concerns the 
long time scale necessary for full strain relaxation. It 
may be possible to modify the time derivative terms in 
the equation of motion to speed up this relaxation. Fur- 
thermore, to make contact with expected results from 
the fracture community (such as the role of the stress 
intensity factor and the fact that crack tips will propa- 
gate without inertia at least until times such that shear 
wave reflections from the boundaries can affect the mo- 
tion) will require much larger systems and much more 
attention to details of the initial conditions. 
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